function[] = Acoustic3rdOctave_v03(first, last, dataset)
global Project;
%------------------------------------------------------------------------------
% 1/3 Octave Band Acoustical Analysis
% This function uses two different methods to calculate the 1/3 octave bands
% Method No. 1: Higher frequencies, direct implementation of filters.
% Method No. 2: Lower frequencies, decimation by series of 3 bands.
%------------------------------------------------------------------------------
%Set variables
pi = 3.14159265358979;
Fref = [20 25 31.5 40 50 63 80 100 125 160 200 250 315 400 500 630 800 1000 1250 1600 2000];
CenterFreq = 1000*((2^(1/3)).^(-17:1:12));  %Exact center freq.
order = 3; %Order of analysis filters.
U = 2^(1/3);
V = 2^(1/6);
Pref = 2.903*10^-9; % Zero dB reference in PSIrms for PSD db calculation.
Fs = Project.Test_Item.Run.Sample_Rate;
w = waitbar(0,'Please wait...');

%------------------------------------------------------------------------------
% For each data channel from 1 to Project.Test_Item.Run.Num_Chan
for ChanNum = 1:Project.Test_Item.Run.Num_Chan
   Project.Test_Item.Run.Channel(ChanNum).Analysis(dataset).Type = 7;
   Project.Test_Item.Run.Channel(ChanNum).Analysis(dataset).Abs = Fref;
   data = Project.Test_Item.Run.Channel(ChanNum).Time_His(first:last);
	DataLength = length(data);

   %Frequency Range
	CenterFreq_high_index = 21;
	CenterFreq_low_index = 1;

    %Check sampling frequencies and issue warnings.
	if (Fs/2) < Fref(CenterFreq_low_index) * 1.5
        error('Sampling Rate too low.')
	elseif (Fs/2) < Fref(CenterFreq_high_index) * 1.5
        disp('Warning: frequency range must be reduced (Sample Rate too Low).');
	    CenterFreq_high_index = find(CenterFreq <= Fs/3,1,'last');
	end % if (Fs/2) < Freq(CenterFreq_low_index) * 1.5 error('Sampling Rate too low.')
	%Compute 'pivot' frequency for multirate filter implementation
	%All filters below Fs/20 will be implemented after a decimation.
	%Fs/20 = 500.

    % CenterFreq_decimate_index: The index of the frequency that marks the lowest frequency
    % where Method No. 1 will be used to calculate the 1/3 octave bands
	CenterFreq_decimate_index = find(CenterFreq <= Fs/20,1,'last');

	OverallRMS = 0;
	%Design filters and compute OverallRMS powers in 1/3-oct. bands.
    %------------------------------------------------------------------------------
    % Method No. 1: Higher frequencies, direct implementation of filters.
    %------------------------------------------------------------------------------
    for index = CenterFreq_high_index:-1:CenterFreq_decimate_index+1
        %Butter is based on a bilinear transformation.
	    f1 = CenterFreq(index) / V;
	    f2 = CenterFreq(index) * V;
        Qr = CenterFreq(index) / (f2-f1);
	    Qd = (pi/2/order)/(sin(pi/2/order))*Qr;
	    alpha = (1 + sqrt(1 + 4*Qd^2))/2/Qd;
	    W1 = CenterFreq(index)/(Fs/2)/alpha;
	    W2 = CenterFreq(index)/(Fs/2)*alpha;
	    [B,A] = butter(order,[W1,W2]);
	    y = filter(B,A,data);

       %Calculate RMS of elements with reference to level Pref.
	    RMS = sqrt(sum(y(1:round(DataLength)).^2)/DataLength);
	    OverallRMS = OverallRMS + RMS;
	    dB = 2*10*log10(RMS/Pref);
	    Project.Test_Item.Run.Channel(ChanNum).Analysis(dataset).Ord(:,index) = dB;
    end % for index = CenterFreq_high_index:-1:CenterFreq_decimate_index+1
    %------------------------------------------------------------------------------
    % End of Method No. 1: Higher frequencies, direct implementation of filters.
    %------------------------------------------------------------------------------

    %------------------------------------------------------------------------------
    % Method No. 2: Lower frequencies, decimation by series of 3 bands.
    %------------------------------------------------------------------------------
	%Finding bandwidth frequencies and coefficients for filters
	NyquistFreq = Fs/2;

	%Calculate coefficients for upper 1/3-oct. band in last octave
	Freq_UpperBand = CenterFreq(CenterFreq_decimate_index);
    f1 = Freq_UpperBand / V;
	f2 = Freq_UpperBand * V;
	Qr = Freq_UpperBand / (f2-f1);
	Qd = (pi/2/order)/(sin(pi/2/order))*Qr;
	alpha = (1 + sqrt(1 + 4*Qd^2))/2/Qd;
	W1 = Freq_UpperBand/(NyquistFreq/2)/alpha;
	W2 = Freq_UpperBand/(NyquistFreq/2)*alpha;
   [Bu,Au] = butter(order,[W1,W2]); % Coefficients for upper 1/3-oct. band

	%Calculate coefficients for center 1/3-oct. band in last octave
	Freq_CenterBand = CenterFreq(CenterFreq_decimate_index)/U;
	f1 = Freq_CenterBand / V;
	f2 = Freq_CenterBand * V;
	Qr = Freq_CenterBand / (f2-f1);
	Qd = (pi/2/order)/(sin(pi/2/order))*Qr;
	alpha = (1 + sqrt(1 + 4*Qd^2))/2/Qd;
	W1 = Freq_CenterBand/(NyquistFreq/2)/alpha;
	W2 = Freq_CenterBand/(NyquistFreq/2)*alpha;
	[Bc,Ac] = butter(order,[W1,W2]); % Coefficients for center 1/3-oct. band

	%Calculate coefficients for lower 1/3-oct. band in last octave
	Freq_LowerBand = CenterFreq(CenterFreq_decimate_index)/(U^2);
	f1 = Freq_LowerBand / V;
	f2 = Freq_LowerBand * V;
	Qr = Freq_LowerBand / (f2-f1);
	Qd = (pi/2/order)/(sin(pi/2/order))*Qr;
	alpha = (1 + sqrt(1 + 4*Qd^2))/2/Qd;
	W1 = Freq_LowerBand/(NyquistFreq/2)/alpha;
	W2 = Freq_LowerBand/(NyquistFreq/2)*alpha;
	[Bl,Al] = butter(order,[W1,W2]); % Coefficients for lower 1/3-oct. band

    index = CenterFreq_decimate_index;
    while index >= CenterFreq_low_index + 2
	    data = decimate(double(data),2);
	    DataLength = length(data);

       y = filter(Bu,Au,data);
       %Calculate RMS of elements with reference to level Pref.
	    RMS = sqrt(sum(y(1:round(DataLength)).^2)/DataLength);
	    dB = 2*10*log10(RMS/Pref);
	    r1 = dB;
	    Project.Test_Item.Run.Channel(ChanNum).Analysis(dataset).Ord(:,index) = dB;

       y = filter(Bc,Ac,data);
       %Calculate RMS of elements with reference to level Pref.
	    RMS = sqrt(sum(y(1:round(DataLength)).^2)/DataLength);
	    dB = 2*10*log10(RMS/Pref);
	    r2 =  dB;
	    Project.Test_Item.Run.Channel(ChanNum).Analysis(dataset).Ord(:,index-1) = dB;

       y = filter(Bl,Al,data);
       %Calculate RMS of elements with reference to level Pref.
	    RMS = sqrt(sum(y(1:round(DataLength)).^2)/DataLength);
	    dB = 2*10*log10(RMS/Pref);
	    r3 = dB;
	    Project.Test_Item.Run.Channel(ChanNum).Analysis(dataset).Ord(:,index-2) = dB;

       index = index - 3;
	    OverallRMS = r1 + r2 + r3 + OverallRMS;
    end % while index >= CenterFreq_low_index + 2;
    %------------------------------------------------------------------------------
	 % End of Method No. 2: Lower frequencies, decimation by series of 3 bands.
    %------------------------------------------------------------------------------

    Project.Test_Item.Run.Channel(ChanNum).Analysis(dataset).FPoint = first;
    Project.Test_Item.Run.Channel(ChanNum).Analysis(dataset).LPoint = last;
    
% 	f =  Project.Test_I7tem.Run.Channel(ChanNum).Analysis(dataset).Abs(:,CenterFreq_low_index:CenterFreq_high_index);
% 	p =  Project.Test_Item.Run.Channel(ChanNum).Analysis(dataset).Ord(:,CenterFreq_low_index:CenterFreq_high_index);

	waitbar(ChanNum/Project.Test_Item.Run.Num_Chan);
end % for ChanNum = 1:Project.Test_Item.Run.Num_Chan
close(w);
